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Theoretical  Studies  of  Oxygen  Reduction  and  Proton  Transfer  in  SOFCs  and  Nerve 

Agents  on  Selected  Surfaces 
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PART  A:  Oxygen  Reduction  in  SOFCs 
I.  Computational  Methods 

All  calculations  were  performed  in  Gaussian09  suite  of  quantum  programs.  Local  minimum 
structures  were  optimized  using  the  generalized  gradient  approximation  (GGA)  based  hybrid 
functional  of  b31yp  combining  with  the  all-electron  basis  sets  and  polarization  functions  of  6- 
31G(d).  The  vibrational  frequencies  were  calculated  at  the  same  theoretical  level,  ensuring  that 
the  reactant  and  product  structures  are  located  at  a  true  minimum  on  the  potential  energy  surface 
(PES)  and  the  transition  state  is  a  first  order  saddle  point  on  PES  as  well.  The  atomic  polar  tensor 
(APT)  charges  were  computed  through  frequency  calculations.  All  resulting  energies  include  the 
zero-point  energy  (ZPE)  corrections.  In  addition,  variable  basis  sets  and  methods  were  applied  to 
verify  the  accuracy  and  reliability  of  the  B3LYP/6-31G(d)  method,  including  6-311G++(d), 
CBSB7,  LANL2DZ,  MP4  and  CCSD(T). 


Gas  C  hannel 


Electrolyte 


Figure  1:  Illustration  of  Oxygen  Reduction  in  SOFCs  (Left:  No  MC,  Right:  With  MC). 


C042  C052 

Figure  2:  Optimized  structures  of  CO42'  and  CO52'  at  the  B3LYP/6-31G  (d)  level. 


II.  Results  and  Discussions 
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2.1  Structure  of  CO42',  CO52'  and  carbonate  clusters 

The  structures  of  CO42"  and  CO52"  were  optimized  at  the  B3LYP/6-31G(d)  level,  as  displayed 
in  Fig.  2.  The  calculated  bond  lengths  and  bond  angles  in  both  CO42'  and  CO52'  are  in  excellent 
agreement  with  the  previous  reported  values.  To  test  the  reliability  of  B3LYP,  we  re-optimized 
the  structure  of  CO52'  using  MP2  with  the  same  basis  set.  All  structural  parameters  are  listed  in 
Table  1.  The  difference  between  the  results  from  two  methods  ranges  0.4-2.9%,  showing  very 
good  agreement.  On  the  other  hand,  the  energy  change  for  C032'+02— >C052'  is  calculated  to  be  - 
105.5  kJ/mol  and  -87.3  kJ/mol  by  B3LYP  and  CCSD(T),  respectively.  Similarly,  the  formation 
energy  of  C042'  (C032'+1/202^C042‘ )  is  -9.8  kJ/mol  and  -5.4  kJ/mol  by  B3LYP  and  CCSD(T), 
respectively.  All  testing  results  have  convinced  the  reliability  of  B3LYP  in  treating  molten 
carbonate  systems. 

Table  1  The  structural  parameters  of  CO52'  obtained  from  B3LYP  and  MP2  (Distance  in 
angstroms,  and  angle  in  degrees). 


01-02 

02-03 

03-C 

01-02-03 

02-03-C 

MP2  1.316 

1.957 

1.380 

110.4 

114.2 

B3LYP  1.321 

1.900 

1.358 

110.8 

116.3 

Diff.%  0.4% 

2.9% 

1.6% 

0.4% 

1.7% 

Fig.  3  shows  the  optimized  structures  of  (Li2C03)4,  (Na2C03)4,  and  (K2C03)4  at  the  b31yp/6- 
31g(d)  level.  As  seen  from  each  picture,  each  alkali  metal  atom  (M)  is  bonded  to  three  carbonate 
oxygen  atoms,  while  each  carbonate  oxygen  atom  is  connected  to  two  M  atoms.  This 
configuration  is  same  as  those  in  their  crystal  structures  of  bulk  Li2C03,  Na2C03,  and  K2C03.  In 
addition,  the  average  bond  length  between  alkali  metal  atom  and  oxygen  (<4i-o)  in  (Li2C03)4, 
(Na2C03)4  and  (K2C03)4  is  calculated  to  be  1.968  A,  2.305  A  and  2.666  A,  respectively,  which 
agrees  well  with  experimental  values  of  1.960  A,  2.377  A  and  2.870  A,  respectively.  Our  recent 
ab  initio  molecular  dynamics  (AIMD)  study  confirms  that  the  volume  expansion  of  Li2C03  is 
only  3%  at  the  temperature  of  1300  K  (Calculated  melting  point  -1000  K),  implying  that  the  Li- 
O  bond  distance  will  not  change  significantly  in  the  molten  lithium  carbonate.  All  evidences 
above  give  us  confidence  in  using  (M2C03)4  to  describe  the  structure  of  molten  carbonates.  More 
importantly,  such  cluster  models  are  more  affordable  when  treating  such  large  molecular  systems. 


4 1 

9*  „  _  * 


(Li2C03)4 
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{1^003)4 


Figure  3:  Optimized  structures  of  (M2C03)4  (M=Li,  Na,  K)  at  the  B3LYP/6-31G  (d)  level. 
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2.2  Oxygen  dissociation  in  the  (Li2C03)4  cluster 

Fig.  4  shows  the  structural  evolution  of  oxygen,  the  beginning,  transitioning  and  ending  states, 
in  the  process  of  dissociation  in  a  lithium  carbonate  cluster  together  with  the  energy  profile.  As 
seen  from  Fig.  4  (a),  the  inserted  figure  is  the  structure  with  minimum  energy  optimized  from  the 
(Li2C03)402  cluster.  It  is  found  that  the  added  oxygen  molecule  bonds  with  one  carbonate  anion 
to  form  CO52",  which  is  circled  by  green  dashed  line.  The  existence  of  CO52'  in  lithium  carbonate 
agrees  with  the  previous  assumption  and  proves  that  oxygen  is  captured  by  carbonate  ions  in  MC. 
From  the  amplified  figure,  the  bond  distance  of  025-026  is  1.427  A,  which  is  0.106  A  longer 
than  that  of  the  free  CO52'  in  the  gas  phase  (Fig.  2).  However,  the  026-010  bond  is  1.474  A  only 
and  0.426  A  shorter  than  in  the  single  CO52".  The  APT  charge  of  025  and  026  is  calculated  to  be 
0.09e  and  -0.69e,  respectively.  The  negative  charge  is  approximately  75%  from  010  and  25% 
from  surrounding  Li  atoms.  Because  026  is  the  end  oxygen  and  connected  to  three  Li  atoms,  it  is 
not  surprising  that  the  electron  will  be  pulled  towards  026  by  Coulomb  force.  The  electron 
density  is  likely  injected  into  the  anti-bonding  orbital  of  025-026,  causing  the  increase  of  the 
bond  distance.  Both  calculated  bond  distances  and  charges  indicate  that  the  O2-CO32'  interaction 
is  of  covalent  bonding  and  the  oxygen  molecule  is  weakened  by  binding  to  carbonate  ion. 
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Figure  4:  Structures  of  beginning,  transitioning  and  ending  state  in  (L^CC^Ch  as  well  as 
PES.  (Distance  in  A;  gray,  red,  and  purple  ball  represents  C,  O,  and  Li,  respectively.) 
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As  025  leaves  026  and  moves  towards  another  carbonate  oxygen  atom  016.  A  new  025-016 
bond  is  forming  while  025-026  is  breaking.  This  leads  to  a  transition  state  (TS)  located  on  the 
PES  (Fig.  3b).  At  TS,  the  distance  between  025  and  026  is  enlarged  to  1.934  A,  while  the 
distance  between  010  and  026  is  0.036  A  shorter.  During  this  process,  the  system  experiences 
an  energy  barrier  of  197.9  kJ/mol.  Finally,  025-026  bond  is  completely  broken  with  a  distance 
of  -3.0  A,  and  025-016  bond  is  formed.  The  PES  calculated  by  intrinsic  reaction  coordinate 
(IRC)  connects  the  starting  structure  with  CO5  ',  the  TS  and  the  final  structure  with  two  CO4  ' 
(Fig.  3c).  This  gives  a  full  picture  of  oxygen  dissociation  in  lithium  carbonate.  Obviously,  with 
the  aid  of  carbonate,  oxygen  dissociation  is  facilitated. 


(c)  Product 


(d)  Energy  Barrier 


Figure  5:  Structures  of  beginning,  transitioning  and  ending  state  in  (Na2C03)402  as  well  as 
PES.  (Distance  in  A;  gray,  red,  and  purple  ball  represents  C,  O,  and  Na,  respectively.) 


2.3  Oxygen  dissociation  in  the  (Na2<303)4  cluster 

Fig.  5  shows  the  structures  of  the  beginning,  transitioning  and  ending  states  for  oxygen 
dissociation  in  a  sodium  carbonate  cluster  as  well  as  the  energy  profile.  Similarly,  the  oxygen 
molecule  of  025-026  was  introduced  to  the  (Na2C03)4  cluster  and  a  bond  was  formed  between 
026  and  04.  The  formed  local  structure  is  also  CO52'  same  as  in  (Li2C03)4  (green  dashed  line  in 


TECH  REPORT  2012-2014 


the  insert  figure).  The  bond  length  of  04-026  and  025-026  is  1.578  A  and  1.391  A,  respectively. 
The  APT  charge  of  025  and  026  is  calculated  to  be  0.2 le  and  -0.60e,  respectively.  As  seen  from 
Fig.  4  (b),  025  moves  away  from  026  towards  06,  initiating  a  dissociation  process.  At  TS,  the 
distance  of  025-026  is  significantly  elongated  to  1.986  A  and  the  04-026  bond  length  is 
slightly  shortened  by  0.114  A.  At  the  end,  025  bonds  with  06  to  form  CO42',  and  the  oxygen 
dissociation  is  completed.  The  calculated  PES  indicates  that  the  oxygen  dissociation  reaction  has 
an  energy  barrier  of  1 16.7  kJ/mol. 


(c)  Product  (d)  Energy  barrier 


Figure  6:  Structures  of  beginning,  transitioning  and  ending  state  in  (K2C03)402  as  well  as 
PES.  (Distance  in  A;  gray,  red,  and  purple  ball  represents  C,  O,  and  K,  respectively.) 


2.4  Oxygen  dissociation  in  the  1X2003)4  cluster 

Fig.  6  shows  the  structural  features  of  the  beginning,  transitioning  and  ending  states  for 
oxygen  dissociation  in  a  potassium  carbonate  cluster  as  well  as  the  energy  profile.  Similarly, 
025-026  is  bonded  to  08  as  a  diatomic  oxygen.  The  08-025  and  025-026  bond  length  is 
calculated  to  be  1.450  A  and  1.408  A,  respectively.  The  APT  charge  on  025  and  026  is  0.1  le 
and  -0.69  e,  respectively.  As  026  moves  away  from  025  towards  to  03,  the  025-026  bond 
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starts  breaking,  while  03-026  bond  starts  forming.  At  TS,  025-026  reaches  a  distance  of  2.000 
A.  The  calculated  energy  barrier  for  the  oxygen  dissociation  is  170.3  kJ/mol.  Finally,  the  oxygen 
dissociation  is  finished  as  the  026-03  bond  is  shortened  to  1.466  A  and  025-026  elongated  to 
3.470  A.  The  reaction  path  on  PES  was  verified  by  IRC  calculations. 

2.5  Mechanism  of  Oxygen  Reduction  by  MC  in  SOFCs 

For  the  reaction  of  oxygen  reduction  (ORR)  in  SOFCs,  the  adsorption  and  dissociation  of 
oxygen  on  the  surface  of  cathode  is  of  great  importance.  Without  MC,  the  reaction  proceeds 
through  superoxide  and  peroxide  intermediates  on  the  surface  of  solid  oxide  cathode  37 .  The 
adsorption  of  oxygen  on  the  ionic  solid  surface  is  limited  to  the  surface  defect  sites  only.  With 
the  presence  of  MC,  the  formation  of  CO52'  was  confirmed  in  three  alkali  molten  carbonates, 
indicating  a  chemisorption  of  gas  oxygen  on  the  surface  of  MC  modified  cathode.  The  binding 
energy  was  estimated  to  be  101.7  kJ/mol.  The  strong  interaction  of  oxygen  with  MC  implies  a 
facilitated  ORR  process  with  the  aid  of  MC  on  the  cathode  of  SOFCs. 

Once  CO52"  is  formed,  the  next  step  is  to  react  with  another  CO32"  to  form  two  CO42'.  This  step 
is  rate-limiting  on  PES  (Fig.  7).  The  calculated  energy  barrier  is  197.9,  116.7,  and  170.3  kJ/mol 
for  lithium,  sodium  and  potassium  molten  salt,  respectively.  The  energy  change  for  this  step  is 
calculated  to  be  +103.3,  +55.6  and  +68.2  kJ/mol,  respectively.  This  fairly  agrees  with  the 
calculated  value  of  +86.2  kJ/mol  for  CO52"  +  CO32'  — »  2CO42"  in  gas  phase  at  the  B3LYP/6- 
31G(d)  level.  Surprisingly,  it  was  found  that  sodium  carbonate  has  the  lowest  activation  energy 
and  also  smallest  energy  change  for  oxygen  dissociation.  This  agrees  with  the  experimental 
observations  that  the  cathodes  modified  by  Li2C03-Na2C03  (52:48  mol  %)  in  SOFCs  have  a 
better  performance  than  those  by  Li2C03-K2C03  (62:38  mol  %)  38.  Table  2  gives  the  bond 
distance  of  C-O  and  M-0  in  the  selected  carbonates,  the  C-0  bond  is  almost  same  in  all  three 
salts,  but  the  M-O  bond  distance  increases  from  Li,  Na  to  K  due  to  the  change  of  atomic  size. 
The  APT  charge  on  Li,  Na  and  K  is  calculated  to  be  0.751  e,  0.769  e  and  0.840  e,  respectively, 
showing  good  consistency  with  the  electronegativity.  The  electrostatic  potential  energy  is 
described  by  QMm-o,  and  the  calculated  value  is  0.382  e/A,  0.334  e/A  and  0.315  e/A, 
respectively.  Lower  potential  energy  and  lighter  the  ions,  the  impendence  for  the  movement  of 
ions  in  MC  will  be  lower  and  so  the  energy  barrier  be  for  oxygen  dissociation  and  diffusion  in 
MC.  This  can  partially  explain  the  lower  energy  barrier  in  sodium  carbonate  than  the  other  two 
salts.  However,  more  details  based  on  MD  simulations  will  be  reported  in  near  future. 


Table  2  Selected  bond  distances  (dc-o  and  dM-o)>  APT  charges  (Qm),  and  electrostatic 
potentials  (QA/m-o)- 


Clusters 

dc-o 

Jm-O 

Qm  QA/m-o 

(Li2C03)402 

1.299 

1.968 

0.751 

0.382 

(Na2C03)402 

1.299 

2.305 

0.769 

0.334 

(K2C03)402 

1.302 

2.666 

0.840 

0.315 

*  Bond  distance  in  A  and  charge  in  e,  all  values  are  average. 
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If  the  reaction  is  approximated  as  Eqn.  5,  then  CO5  2‘  can  be  treated  as  an  intermediate  and  the 
steady-state  theory  can  be  applied.  The  reaction  rate  is  described  by  Eqn.  6.  The  effective 
activation  energy  for  the  overall  reaction  is  estimated  to  be  96.2,  15.1,  and  68.6  kJ/mol  in  lithium, 
sodium  and  potassium  molten  carbonate,  respectively.  The  pseudo  one-step  reaction  is 
exothermic  and  favored  by  chemical  thermodynamics  for  Na  and  K,  but  slightly  endothermic  for 
Li.  From  the  view  of  chemical  kinetics,  lithium  carbonate  is  not  a  favorable  molten  salt  for  ORR, 
but  lithium  carbonate  can  significantly  reduce  the  melting  point  of  sodium  and  potassium 
carbonates  when  eutectic  ratios  are  adopted.  In  Eqn.  6,  the  concentration  of  CO32'  is  a  constant, 
therefore  the  reaction  rate  is  solely  determined  by  the  pressure  of  oxygen.  Increasing  the  oxygen 
pressure  will  improve  the  kinetics  of  the  cathode  in  SOFCs. 

2C02-  +  02  <-^->C02~  +  CO2'  — — — »2C02'  (5) 

k- 1 

Reaction  Rate  =  [02 ]  [CO2  ]2  (6) 

Given  that  01-02-03  in  CO52'  (Fig.  1)  is  ozone-like  and  possibly  has  considerable  biradical 
character,  it  is  necessary  to  verify  the  reliability  of  B3LYP  in  this  study.  Previous  studies  by 
Gerber  and  co-workers  on  Criegee  Intermediates  (CIs)  have  proven  that  MP4  and  CCSD(T)  are 
reasonably  reliable  in  treating  such  molecular  systems.  Therefore,  we  have  calculated  single 
point  energies  using  variable  basis  sets  and  methods.  Calculated  energy  barriers  for  oxygen 
dissociation  are  summarized  in  Table  3.  Only  small  changes  occur  when  the  basis  set  is  changed 
from  6-31G(d),  6-311G++(d)  to  CBSB7,  all  using  B3LYP.  Considering  the  computing  cost  with 
MP4  and  CCSD(T),  we  had  to  decrease  the  basis  set  to  LANL2DZ.  The  use  of  effective  core 
potentials  (ECPs)  has  largely  reduced  the  computing  time.  The  calculated  energy  barriers  by 
MP4  are  higher  than  those  by  B3LYP,  about  10%  in  (K2C03)4  and  17%  in  (Na2C03)4,  but  only  1% 
in  (Li2C03)4.  However,  an  excellent  agreement  was  observed  between  B3LYP  and  CCSD(T). 
The  comparative  results  have  verified  that  the  results  by  the  B3LYP/6-31G(d)  method  for  the 
oxygen  dissociation  in  (M2C03)4  (M=Li,  Na,  K)  clusters  are  accurate  and  the  method  is  reliable 
for  similar  studies  in  future. 


Figure  7:  PES  of  Oxygen  Dissociation  in  the  (M2C03)4  clusters,  M=Li,  Na,  and  K. 
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Table  3  The  activation  energy  (kj/mol)  of  oxygen  dissociation  in  (M2C03)4  (M=Li,  Na,  K). 


B3LYP/6-31G(d)  B3LYP/6-311++G(d)  B3LYP/CBSB7 


(Li2C03)4 

204.2/1 97.9(ZPE) 

198.7 

197.5 

(Na2C03)4 

1 15.9/1 16.7(ZPE) 

115.9 

115.1 

(K2co3)4 

171.1/170.3(ZPE) 

159.8 

159.4 

B3LYP/LANL2DZ 

CCSD(T)/LANL2DZ 

MP4/LANL2DZ 

(Li2C03)4 

197.1 

199.6 

194.6 

(Na2C03)4 

87.9 

103.3 

95.4 

(K2co3)4 

141.0 

154.0 

146.0 

III.  Summary  of  Part  A 

We  have  reported  a  DFT  study  of  oxygen  dissociation  process  in  selected  molten  carbonates. 
It  was  found  that  the  adsorption  of  oxygen  to  molten  carbonate  is  of  chemical  type,  and  leads  to 
the  formation  of  CO52'  in  MC,  which  was  confirmed  for  the  first  time  by  DFT  calculations.  The 
energy  barrier  for  its  dissociation  is  calculated  to  be  197.9,  116.7,  and  170.3  kJ/mol  at  the 
B3LYP/6-31G(d)  level  in  the  (M2C03)4  cluster,  M  =  Li,  Na,  and  K,  respectively.  If  the  reaction 
of  02  +  2CO32"  — »  2C042~  is  approximated  as  a  one-step  reaction,  the  activation  energy  is 
estimated  to  be  96.2,  15.1,  and  68.6  kJ/mol,  respectively.  The  reaction  rate  is  the  first  order  to  the 
pressure  of  oxygen.  To  our  surprise,  sodium  carbonate  has  the  lowest  energy  barrier  for  oxygen 
dissociation  in  molten  sodium  carbonate,  which  is  consistent  to  the  recent  experimental  findings. 
It  is  very  clear  that  the  molten  carbonate  salt  has  directly  participated  in  the  ORR  process  and 
plays  an  important  role  as  a  catalyst.  The  oxygen  reduction  has  been  facilitated  by  MC  and 
enhanced  cell  performance  has  been  observed.  In  addition,  the  B3LYP/6-31G(d)  method  is 
verified  to  be  accurate  and  reliable  in  treating  such  molten  carbonate  salt  systems. 
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Theoretical  Studies  of  Oxygen  Reduction  and  Proton  Transfer  in  SOFCs  and  Nerve 

Agents  on  Selected  Surfaces 

Changyong  Qin,  Benedict  College,  Columbia,  SC  29204  (PI) 


PART  B:  Mechanism  of  Proton  Transfer  in  SOFCs 
I.  Computational  Methods 

1.1  Static  DFT  Calculations  in  Gaussian  09 

Static  DFT  calculations  were  carried  out  at  the  B3LYP/6-31G(d,p)  level  as  implemented  in  the 
Gaussian  09  package.  Geometries  of  ground-state  and  transition  state  (TS)  of  [(Li2C03)8H]+ 
clusters  were  fully  optimized.  To  confirm  each  optimized  stationary  point,  the  vibrational 
frequencies  were  examined  at  the  same  theoretical  level,  ensuring  that  the  reactant  and  product 
structures  are  local  minimum  while  the  transition  state  is  a  first  order  saddle  point  on  the 
potential  energy  surface  (PES).  To  quickly  and  correctly  search  out  the  TS,  we  first  performed 
the  SCAN  calculations  and  then  performed  TS  calculations  using  the  geometry  with  the  highest 
energy  .  Finally,  the  intrinsic  reaction  coordinate  (IRC)  calculations  were  carried  out  to  judge  the 
reaction  direction. 

1.2  Dynamic  DFT  Calculations  in  VASP 

Using  VASP  5.3  package,  first-principles  molecular  dynamics  (FPMD)  calculations  were 
performed.  The  structure  of  [(Li2C03)8H]+  cluster  optimized  at  the  B3LYP/6-31G(d,p)  level  was 
used  as  the  initial  structure  of  FPMD  simulations.  All  FPMD  simulations  were  conducted  with 
the  NVT  canonical  ensemble.  PAW-PBE  potentials  supplied  with  the  VASP  package  were  used 
for  H  (ultrasoft),  Li  (slpO),  C  (s2p2),  and  O  (s2p4).  The  Verlet  algorithm  was  integrated  with 
Newton’s  equations  of  motion  at  a  time  step  of  2  fs  for  a  total  simulation  time  of  6  ps,  i.e.,  3000 
steps.  The  frequency  of  the  temperature  oscillations  was  controlled  by  the  Nose  mass  during  the 
simulations.  Energy  cutoff  of  500  eV  and  a  lxlxl  k-point  mesh  were  used. 

For  the  simulations  of  proton  transfer  in  the  bulk  lithium  molten  carbonate,  the  supercell  of 
1x2x2  was  first  melted  by  virtually  heated  to  3000  K  within  the  NVT  ensemble.  Then,  a  melt 
configuration  was  randomly  chosen  to  test  the  optimal  lattice  constant  while  considering  the 
change  in  system  volume  induced  by  the  phase  transition.  Then,  the  melt  configuration  was 
subjected  to  a  cooling  process  at  a  rate  of  2.67xl013  K/s  to  obtain  a  molten  carbonate  structure  at 
1100  K,  1200  K,  and  1300  K.  Finally,  we  randomly  put  the  proton  in  the  lithium  molten 
carbonate  and  simulated  for  another  6  ps  at  the  temperatures  of  1100  K,  1200  K,  and  1300  K, 
respectively. 

II.  Results  and  Discussion 

2.1  Intramolecular  Transfer  of  Proton  in  the  Carbonate  Ion 

First,  we  studied  the  proton  migration  in  intra-(C03)2'  ion.  As  shown  in  Fig.l,  during  the  transfer, 
the  bond  between  Ol  and  H+  is  first  broken,  followed  by  H+  moving  toward  the  mirror  position 
between  Ol  and  02,  namely  the  TS.  The  bond  between  02  and  H+  is  then  formed,  completing 
the  H+-transfer.  We  find  that  the  energy  barrier  for  H+-transfer  between  Ol  and  02  is  22.3 
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kcal/mol.  In  addition,  to  test  the  reliability  of  B3LYP,  the  CCSD(T)  method  with  the  same  basis 
set  was  used  to  compute  the  single  point  energies.  The  activation  energy  computed  at 
CCSD(T)/6-31G  (d,p)  level  is  25.5  kcal/mol,  which  is  only  slightly  higher  than  that  obtained 
from  B3LYP/6-31G  (d,p),  about  3.2  kcal/mol.  The  comparative  activation  energy  has  convinced 
the  B3LYP  method  is  reliable  in  treating  such  systems. 


- B3LYP 

- CCSD(T) 


O _ r,- 

Reactant 


25.5  kcal/mol 


Product 


Fig.  1  The  structures  of  reactant,  transition  state  (TS),  and  product  as  well  as  the  relative 
energy  for  proton  transfer  in  intra-(C03)2'  ion.  The  proton,  carbon,  and  oxygen  atoms  are 
shown  as  white,  gray,  and  red  balls,  respectively. 

Second,  we  studied  the  proton  transfer  in  intra-carbonate  and  inter-carbonate  in  (Li2C03)8  cluster. 
Fig.  2(a)  and  Fig.  2(b)  indicates  the  proton  transfer  in  intra-carbonate  on  the  surface  of  (LioCCfdg 
cluster  and  inside  the  (Li2C03)g  cluster,  respectively.  As  seen  from  Fig.  2,  both  migrations 
experience  the  H+-0  bond  broken  and  reforming,  and  the  energy  barriers  are  46.8  kcal/mol  and 
49.0  kcal/mol,  respectively. 
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(a) 

TS 

46.8  kcal/mol 

\  40.4  kcal/mol 

reactant 

product 

(b) 

TS 

49.0  kcal/mol 

\  41 .0  kcal/mol 

reactant 

product 

Fig.  2  The  structures  of  reactant,  transition  state  (TS),  and  product  as  well  as  the  relative 
energy  for  proton  migration  in  intra-carbonate  in  the  (L^CC^s  cluster.  The  proton, 
carbon,  oxygen  and  lithium  atoms  are  shown  as  yellow,  gray,  red,  and  purple  respectively. 

2.1  Intermolecular  Transfer  of  Proton  between  Carbonate  Ions 

Fig.  3(a)  and  Fig.  3(b)  indicates  the  proton  transfer  in  inter-carbonate  on  the  surface  of  (L^CCF)* 
cluster  and  inside  the  (L^CCFlg  cluster,  respectively.  As  seen  from  Fig.  3,  unlike  proton  transfer 
in  intra-carbonate  as  shown  in  Fig.2,  although  migrations  of  proton  transfer  in  inter-carbonate 
experience  the  H+-0  bond  broken  and  reforming,  but  the  energy  barriers  are  decreased  to  1.1 
kcal/mol  and  7.7  kcal/mol,  respectively.  This  is  not  surprising  because  when  the  proton  transfers 
between  carbonate  ions  (inter-carbonate),  each  carbonate  ion  can  adjust  its  position  accordingly, 
decreasing  the  energy  barrier  of  H+-transfer. 


reactant  fS  product 
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TS 

0.8  kcal/mol 

product 

(b) 

7.7  kcal/mol 

reactant 

(a) 

TS  0.02  kcal/mol 

product 

1.1  kcal/mol 

reactant 

Fig.  3  The  structures  of  reactant,  transition  state  (TS),  and  product  as  well  as  the  relative 
energy  for  proton  migration  in  inter-carbonate  in  the  (LiiCO^s  cluster.  The  proton,  carbon, 
oxygen  and  lithium  atoms  are  shown  as  yellow,  gray,  red,  and  purple  respectively. 

Apparently,  comparison  of  the  above  energy  barrier  suggests  that  there  are  no  energetics 
difference  for  proton  transfer  on  the  surface  of  (L^CC^jg  cluster  and  inside  the  (Li2C03)8  cluster, 
and  the  proton  transfer  prefers  inter-carbonate  to  intra-carbonate  in  the  (L^CCfdg  cluster 

2.3  Dynamic  DFT  on  Proton  Transfer  in  (LTCCXdg 

To  simulate  the  proton  migration  features  in  real  molten  carbonate  at  finite  temperatures,  we 
carried  out  the  first-principles  molecular  dynamics  (FPMD)  calculations.  Fig.  4  shows  the 
trajectories  of  proton,  selected  carbon  and  oxygen  atom  in  the  (L^CC^jg  cluster  simulated  at 
1100  K,  1200  K,  and  1300  K  for  6  ps.  We  can  see  that  the  trajectory  of  proton  is  very  dispersed 
while  the  trajectory  of  carbon  is  very  localized,  which  indicates  that  the  proton  transfer  in  the 
(Li2C03)8  cluster  is  quite  fast  with  large  displacement,  but  the  carbon  atom  only  vibrates  at  its 
original  site  and  unsuccessful  jumps  with  much  limited  displacement.  It  should  be  noted  that  the 
trajectory  of  selected  oxygen  atom  is  dispersed  somewhat,  indicating  that  the  oxygen  atom  adjust 
its  position  by  rotation  or  rolling  during  the  proton  transfer  in  inter-carbonate  in  the  (L^CC^jg 
cluster. 

The  mean  square  displacement  (MSD)  is  determined  by  the  ensemble  average: 

M  SD  =<  >=  1  -  £  f>.  (» +  (0, )  -  r„  (1), 

N  ni=0m=0 

where  N  is  the  number  of  ions  in  the  system,  n  is  the  number  of  time  origins,  t  is  the  time,  and  to; 
is  the  initial  timestep  originating  at  time  i.  In  this  work,  at  least  three  FPMD  running  were 
conducted  for  the  same  structure  at  each  temperature,  aiming  to  obtain  an  averaged  MSD  for 
accuracy.  The  diffusion  coefficient  is  calculated  by  fitting  the  slope  of  the  MSD  vs.  time.  To 
obtain  a  rational  diffusion  coefficient,  only  the  second  half  MSD  data  are  included  in  the  linear 
fitting.  Fig.  5  shows  the  calculated  mean  square  displacement  for  proton  transfer  in  the  (L^CCXjg 
cluster  simulated  at  the  temperature  of  1100  K,  1200  K,  and  1300  K  for  6  ps.  Evidently,  the 
MSD  of  proton  increases  linearly  with  time,  indicating  a  fast  proton  migration.  Moreover,  the 
slope  of  MSD  increases  as  the  temperature  increases,  indicating  the  ability  of  proton  diffusion 
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increases  with  the  temperature  increasing. 


Fig.  4  The  proton,  carbon,  and  oxygen  atom  trajectories  in  the  (LhCC^js  cluster  simulated 

at  1100  K,  1200  K,  and  1300  K  for  6  ps. 


Fig.  5  The  average  mean  square  displacement  (MSD)  for  proton  transfer  in  the  (LiiCC^is 
cluster  simulated  at  1100  K,  1200  K,  and  1300  K  for  6ps. 
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Also,  we  can  give  the  Arrhenius  plot  via  Arrhenius  relationship: 

A  F 

D(T)  =  D0(T)exp(-— — )  (4). 

kBT 

Herein,  D0  and  A E  are  the  pre-exponential  and  the  diffusion  barrier,  respectively.  kB  and  T  are 

the  Boltzmann  constant  and  absolute  temperature,  respectively.  By  fitting  to  the  Arrhenius 
relationship  over  the  temperature  range,  the  pre-exponential  and  the  effective  activation  barrier 
can  be  obtained.  Fig.  6  shows  the  Arrhenius  plot  for  proton  transfer  in  the  (Li2C03)s  cluster. 
From  the  slope  of  the  straight  line,  the  activation  energy  is  calculated  to  be  6.4  kcal/mol,  which  is 
in  good  agreement  with  the  energy  barrier  of  proton  transfer  in  inter-carbonate  in  the  (L^COsjg 
clusters.  This  again  supports  that  the  proton  fast  transfer  prefer  in  inter-carbonate  and  suggests 
that  our  cluster  model  is  suitable  to  describe  the  state  of  molten  carbonate. 


103/T(K_1) 

Fig.  6  Arrhenius  plots  for  proton  transfer  in  the  (Li2C03)g  cluster. 

2.4  Dynamic  DFT  on  Proton  Transfer  in  Bulk  Lithium  Carbonate 

Additionally,  we  also  simulate  the  proton  migration  in  bulk  lithium  carbonate  at  finite 
temperatures  by  using  FPMD.  Since  the  system  volume  is  adjusted  due  to  the  phase  transition 
from  crystal  phase  to  molten  phase,  and  only  constant  volume  MD  are  possible  at  the  moment 
for  the  standard  first-principles  molecular  dynamics  in  the  VASP  program.  So,  we  need  to  test 
the  lattice  constant  for  molten  Li2C03  structure.  As  described  in  theoretical  details,  we  randomly 
chose  a  melt  configuration,  and  then  relax  it  for  20  ps  (NVT)  at  different  volume.  Fig.  7  shows 
the  average  free  energy  versus  the  different  lattice  constant  scale  factors.  We  can  see  that  there  is 
small  influence  of  volume  on  the  free  energy.  The  best  lattice  constant  scale  factor  is  1.010, 
which  means  the  volume  expands  3.0  %  after  melting.  In  our  following  FPMD  simulations,  the 
lattice  parameters  of  all  structures  multiply  by  1.010  to  obtain  the  accurate  potential  energy  and 
free  energy.  Fig.  8  displays  the  structures  of  bulk  lithium  carbonate  before  and  after  melting.  The 
molten  lithium  carbonate  is  simulated  at  3000  K  for  20  ps.  Obviously  the  structures  of  bulk 
lithium  carbonate  after  melting  can  represent  the  state  of  molten  lithium  carbonate. 
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Lattice  constant  seal  factors 

Fig.  7  The  average  free  energy  versus  the  lattice  constant  scale  factors. 


Fig.8  The  structures  of  crystal  (a)  and  molten  phase  (b),  the  latter  is  simulated  at  3000  K. 

Then,  we  randomly  put  one  proton  into  the  molten  lithium  carbonate  and  simulated  at  1100  K, 
1200  K,  and  1300  K  for  6ps,  respectively.  Fig.  9  shows  a  series  of  the  snapshots  of  proton 
positions  at  1200  K  in  a  time  interval  of  1  ps  for  a  total  of  5  ps,  the  white  ball  represents  one 
mobile  proton.  Comparison  of  proton  positions  at  each  snapshot  indicates  a  fast  proton  diffusion 
as  evidenced  by  the  large  change  in  displacement  between  each  snapshot.  The  corresponding 
trajectory  traces  of  proton  as  well  as  C  and  O  ions  are  further  depicted  in  Fig.  10.  It  is  evident  that 
the  proton  transfer  features  like  those  in  the  (LiiCC^jg  cluster:  the  trajectory  of  proton  is  very 
dispersed  while  the  trajectory  of  carbon  is  very  localized,  and  the  trajectory  of  oxygen  is 
somewhat  localized.  This  indicates  that  the  mechanism  of  proton  transfer  in  the  lithium  molten 
carbonate  is  the  same  as  in  the  (L^CC^jg  cluster:  the  proton  fast  transfers  between  the  carbonate 
ions,  at  the  same  time,  the  carbonate  ion  adjust  its  position  accordingly  to  reduce  the  energy 
barrier. 
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Fig.  9  The  snapshots  of  proton  diffusion  in  the  lithium  molten  carbonate  simulated  at  1200 
K.  The  white  ball  indicates  the  proton  (denoted  by  arrow). 


Fig.  10  The  trajectories  of  proton,  carbonate,  and  oxygen  atoms  in  the  lithium  molten 
carbonate  simulated  at  1200  K  for  6  ps. 

Likewise,  we  also  computed  the  mean  square  displacement  (MSD)  and  Arrhenius  plot.  As  shown 
in  Fig.  11  and  Fig.  12,  the  features  of  MSD  are  similar  to  those  in  the  (L^CCDg  cluster:  the  MSD 
of  proton  increases  linearly  with  time  and  the  slope  of  MSD  increases  as  the  temperature 
increases,  indicating  a  fast  proton  migration  in  the  lithium  molten  carbonate.  The  activation 
energy  is  calculated  to  be  2.2  kcal/mol,  which  agrees  well  with  that  in  the  (Li2C03)g  cluster. 


TECH  REPORT  2012-2014 


Such  small  activation  energy  again  indicates  that  the  proton  transfer  in  the  lithium  molten 
carbonate  is  very  facile  and  fast. 


Fig.  11  The  average  mean  square  displacement  (MSD)  for  proton  transfer  in  the  lithium 
molten  carbonate  simulated  at  1100  K,  1200  K,  and  1300  K  for  6ps. 


Fig.  12  Arrhenius  plots  for  proton  transfer  in  the  lithium  molten  carbonate. 

III.  Conclusion  and  Future  Work 

In  summary,  the  proton  transfer  in  molten  lithium  carbonate  has  been  systematically  examined 
using  DFT  calculations.  The  transfer  mechanism  involves  a  concerted  bond  forming  and 
breaking  process.  The  inter-carbonate  pathway  is  preferred  than  the  intra  path  with  calculated 
energy  barrier  of  1-8  kcal/mol  and  40-50  kcal/mol,  respectively.  The  DFT  dynamic  calculations 
confirmed  the  vitality  of  the  inter-carbonate  path.  Calculated  activation  energy  through 
Arrhenius  plotting  is  2.2  and  6.4  kcal/mol  in  MC  cluater  or  bulk  MC  system,  respectively.  The 
analysis  of  atomic  trajectory  shows  high  mobility  of  proton,  while  carbonate  atoms  are  largely 
limited  and  oxygen  atoms  with  possible  rotations  to  accept  the  proton.  Our  next  goal  is  to 
examine  the  proton  transfer  at  the  surface  of  BYZ  and  the  MC/BZY  interface. 
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Theoretical  Studies  of  Oxygen  Reduction  and  Proton  Transfer  in  SOFCs  and  Nerve 

Agents  on  Selected  Surfaces 

Jerry  Whitten,  North  Carolina  State  University,  Raleigh,  NC  27695  (Co-PI) 


PART  C:  Adsorption  of  Nerve  Agents  on  Selected  Surfaces 
I.  Overview  of  Research 

In  this  work,  we  examine  fundamental  processes  involving  the  interaction  of  nerve  agents 
with  solid  surfaces  using  accurate  electronic  structure  methods.  Adsorption  on  several  types  of 
surfaces  found  in  structural  materials  are  considered  for  two  types  of  nerve  agents:  sarin,  a  nerve 
gas,  (see  Fig.  1)  and  VX  type  nerve  agents  .  Although  both  nerve  agents  are  organophosphates, 
they  differ  in  the  ester  linkage  and  in  substituent  groups.  These  differences  are  known  to  lead  to 
different  persistence  times  on  surfaces. 
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Fig.  1.  Sarin  nerve  gas.  An  organo- 
phosphate  nerve  agent  with  an  ester 
linkage  to  an  aliphatic  group. 


The  objective  of  the  work  is  to  describe  surface  interactions  and  low  energy  electronic 
states  at  high  aaccuracy  using  theory  in  order  to  understand  the  details  of  bonding  to  surfaces. 
The  plan  is  to  carry  out  the  study  so  that  factors  that  affect  desorption  energies,  desorption 
kinetics  and  the  spectral  features  of  adsorbed  nerve  agents  can  be  quantified.  Understanding 
exactly  how  these  molecules  interact  with  a  variety  of  solid  surfaces  is  the  primary  goal  of  this 
work. 


The  project  is  divided  into  three  phases: 

1)  Adsorption  of  sarin  on  hydrocarbon  and  oxide  surfaces 

2)  Solvation  of  sarin  adsorbed  on  surfaces  by  water 
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3)  Study  of  the  nerve  agent  VX  and  examination  of  the  low  lying  excited  states  of  sarin 
adsorbed  on  surfaces. 

Phase  1  has  been  completed  and  Phase  2  is  nearing  completion.  The  project  benefitted 
from  the  participation  of  two  sabbatical  leave  visitors  to  the  Whitten  group  from  Athens,  Greece, 
G.  Petsalakis  and  I.  Petsalakis.  The  visitors  had  prior  experience  using  Density  Functional 
Theory  to  describe  molecule-surface  interactions  and  this  permitted  a  comparison  of  DFT  and 
our  many-electron  Cl  calculations  for  sarin  bound  to  hydrocarbon  and  oxide  surfaces. 

Appendix  I  contains  the  manuscript  describing  the  results  from  Phase  1  of  the  project. 
This  manuscript  was  submitted  to  J.  Phys.  Chem.  and  is  currently  undergoing  minor  revision  to 
clarify  questions  raised  in  the  review. 

In  the  Phase  2  solvation  studies,  the  surfaces  are  either  initially  coated  by  water,  see  for 
example  the  snapshot  in  Fig.  2,  followed  by  introduction  of  sarin  at  a  series  of  distances  from  the 
surface,  or  the  sarin  molecule  is  solvated  first,  see  Fig.  3,  and  then  allowed  to  approach  the 
surface.  In  both  cases,  sarin  and  the  surface  compete  for  the  water  molecules  and  the  molecules 
move  optimally  to  solvate  both  sarin  and  the  surface  as  the  molecule  approaches  the  surface.  We 
follow  the  energetics  of  the  system  as  a  function  of  the  distance  of  sarin  from  the  surface 
allowing  the  water  molecules  to  move  to  their  minimum  energy  location  for  each  choice  of  sarin 
position.  The  results  should  lead  to  an  understanding  of  the  effect  of  water  solvation  on  the 
ability  of  a  surface  to  retain  sarin. 


Fig.  2.  Binding  directly  to  water  layer:  16  kcal/mol 
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Fig.  3.  Solvation  of  sarin  near  an  aliphatic  hydrocarbon  surface  by  water. 

II.  Theoretical  Methods 

Two  types  of  electronic  structure  approaches  are  employed:  the  first  one  involves 
configuration  interaction  along  with  a  simplex  optimization  method  for  the  geometry 
optimizations.  The  second  approach  involves  density  functional  theory  (DFT)  calculations  for 
the  ground  electronic  state  and  time  dependent  DFT  (TDDFT)  calculations  for  excited  electronic 
states,  with  details  given  below. 

2.1  Many-electron  Cl  Theory 

The  adsorbate-surface  systems  will  be  described  by  a  configuration  interaction 
embedding  method  that  is  designed  to  give  an  accurate  many-electron  description  of  the 
adsorbate- surface  region  for  both  ground  and  excited  electronic  states.  Many  applications  to 
adsorbates  on  metals  and  oxides  have  been  reported  previously  and  details  can  be  found  in  Ref. 
18-21.  In  the  case  of  small  particles,  there  is  no  need  to  make  approximations  in  the  coupling  of 
the  local  region  to  the  bulk  and  in  this  case  “embedding”  refers  to  the  creation  of  an  electronic 
subspace  that  is  treated  by  configuration  interaction.  A  brief  summary  of  the  theory  is  given 
below. 


Calculations  are  carried  out  for  the  full  electrostatic  Hamiltonian  of  the  system 

N  Q  _  v  N 


N  U  _  y  N  i 
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Wavefunctions  are  constructed  by  self-consistent-field  (SCF)  and  multi-reference  configuration 

V  =  (n!  f  del  (Zf  Z*  ■  ■  ■  zt )  = 


interaction  (Cl)  expansions, 
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Single  and  double  excitations  from  an  initial  representation  of  the  state  of  interest, 
<£>r,  are  carried  out  to  create  a  small  Cl  expansion, 

V1V  =  Or  +  Zijki  A,jjk|  r^ki  Or  =  Om 

Configurations,  ®k,  are  retained  if  the  interaction  with  ®r  satisfies  a  relatively  large  threshold  condition 
|<®k|  H  |  ®r>|  2/  |  Ek  -  Er  |  >  10  '4  a.u. 

Next,  we  refine  the  description  of  the  state  by  generating  a  large  Cl  expansion,^,  by  single  and  double 
excitations  from  all  important  members  of  TV  (those  with  coefficient  >  0.05)  to  obtain 

Tr  =VPr'  +  Zm[Zik  ^ikm  rm®m  +  2ijk|  ^ijk|m  rij4k|®m] 

where  ®m  is  an  important  member  of  TV.  The  additional  configurations  are  generated  by  identifying 
and  retaining  all  configurations,  ®k,  that  interact  with  4V  such  that 

|<®k|  H  |  'IV  >  |  2  /  I  Ek  -  Er  |  >  3x10  7  a.u. 

For  this  small  threshold,  typically  4xl04  configurations  occur  in  the  final  Cl  expansion,  and  the 
expansion  can  contain  single  through  quadruple  excitations  from  an  initial  representation  of  the 
state  Dr.  Contributions  of  configurations  not  explicitly  retained  are  estimated  using  perturbation 
theory. 

In  order  to  accelerate  the  convergence  of  the  Cl  expansion  results  with  respect  to  active 
space  size  for  the  excited  state  calculations,  the  molecular  orbitals  of  the  lowest  triplet  state  for 
each  system  are  first  put  through  a  unitary  transformation.  This  transformation  creates  a  set  of 
occupied  and  virtual  orbitals  which  have  maximal  exchange  interaction  with  the  occupied 
orbitals  within  the  active  space.  The  primary  function  of  this  is  to  allow  accurate  description  of 
the  states  of  interest  in  the  presence  of  less  accurately  described  lower-lying  states.  The  basis  of 
this  argument  lies  in  the  way  other  configurations  are  generated  from  the  states  of  interest.  For  a 
given  excited  state,  Dr,  allowing  all  single  and  double  excitations  from  its  important 
configurations  also  generates  all  lower  lying  states  that  have  significant  interaction  with  Dr. 

Since  these  configurations  are  in  principle  included  in  the  diagonalization,  it  follows  that  the 
excited  state,  Dr,  is  orthogonal  to  and  non-interacting  with  lower  energy  states.  Thus  the 
variational  theorem  for  excited  states  is  satisfied. 

For  the  Cl  computations,  the  interaction  of  sarin  is  described  using  fixed-geometry 
surfaces  sufficiently  large  enough  to  ensure  that  sarin  would  remain  fully  over  the  surface  and 
not  interact  appreciably  with  boundary  atoms.  These  surface  models,  depicted  in  a  subsequent 
figure,  correspond  to  C22H12  for  graphene,  C22H34  for  graphane,  and  Ca^Cfis  for  calcium  oxide. 

Compact,  but  still  flexible  bases  were  used  in  this  study,  and  can  roughly  be  characterized 
as  double-zeta  plus  polarization.22  These  are  described  as  (10s5pld)— >[4s2pld]  for  C  and  O, 
(12s8pld)— >[6s4pld]  for  P,  (10s5p)— >[4s2p]  for  F,  and  (5slp)— >[2slp]  for  H.  For  the  graphane 
surface  only  the  pz  orbital  was  added  to  each  H  atom  to  ensure  sufficient  bond  polarization. 
Graphene  hydrogen  atoms  did  not  use  any  p  functions,  while  the  carbon  atoms  had  additional 
diffuse  s  and  pz  functions  added  with  exponents  of  0.08.  For  calcium  oxide  the  surface  should 
show  little  reaction  to  the  addition  of  sarin  and  so  a  (22sl2p)— >[5s2p]  basis  for  Ca  and  a 
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(12s5p)— >[2slp]  basis  for  O  were  used.  This  basis  was  generated  by  minimizing  the  total  energy 
of  CaiOzt,  and  is  given  in  the  supplemental  information. 

2.2  Simplex  Optimization  Method 

The  Cl  geometry  optimizations  are  carried  out  by  the  Nelder-Mead  simplex  procedure.23 
The  optimization  allows  variation  of  the  geometry  of  adsorbates  on  the  surfaces  of  interest.  In 
the  simplex  procedure,  the  parameters  that  define  the  geometry  are  taken  as  components  of  a  so- 
called  vertex  or  n-tuple.  The  energy  is  calculated  for  that  vertex  and  a  systematic  variation  of  the 
parameters  is  carried  out  in  a  way  that  eliminates  the  worst  vertex  after  each  iteration.  The 
procedure  is  widely  used  in  applied  mathematics  and  is  surprisingly  robust.  It  has  an  advantage 
in  that  energy  gradients  are  not  required  and  the  energy  minimization  can  be  carried  out  directly 
at  the  Cl  level  of  theory  since  only  a  computer  code  capable  of  computing  the  energy  is  needed. 
The  latter  point  is  important  when  electron  correlation  has  a  significant  effect  on  the  geometry, 
e.g.,  stretched  bonds  or  pre-dissociation.  Geometry  optimizations  at  the  Cl  level  allowed  all 
degrees  of  freedom  of  sarin  to  vary,  while  keeping  the  model  surfaces  locked  at  a  standard 
geometry. 

2.3  DFT  and  TDDFT  Calculations 

Complementary  to  the  Cl  calculations  described  in  Section  2.1  DFT24  and  TDDFT25 
calculations  employing  the  M062X  functional  ’  and  two  types  of  basis  sets,  6-3  lG(d,p)  and  6- 
311G+(d,p)  (i.e.  including  diffuse  functions)  have  been  carried  out  on  sarin  adsorbed  on 
graphene,  graphane  and  CaO,  as  in  the  Cl  calculations  but  using  slightly  different  models  for  the 
surfaces.  The  smaller  basis  set  was  employed  for  larger  models  of  the  surfaces,  cf.,  structures  (a) 
and  (d)  for  graphene  (C66H22)  and  graphane  (C66H88),  respectively,  and  structure  (g)  for  CaO  in 
Figure  3.  The  larger  basis  set  was  employed  for  smaller  models  of  the  graphene  (C28H14)  and 
graphane  (C28H42)  surfaces,  (b)  and  (e),  respectively,  in  Figure  3,  and  also  for  the  CaO  (Cai 50 1 5) 
surface,  (h).  Binding  energies  of  sarin  on  the  three  surfaces  have  been  determined,  taking  into 
account  the  basis  superposition  error.  In  addition,  excited  states  and  absorption  spectra  of  the 
adsorption  complexes  of  sarin  on  graphene  and  graphane  have  been  calculated  and  compared  to 
the  absorption  spectrum  of  free  sarin  obtained  by  TDDFT/M062X/6-31 1+G(d,p).  For  further 
comparison,  the  absorption  spectrum  of  isolated  sarin  has  been  calculated  by  TDDFT  employing 
nine  additional  functionals  (specified  in  Figure  4)  and  the  6-311+G(d,p)  basis  set.  Corresponding 
data  has  been  generated  by  Cl  calculations  for  free  sarin  as  well  as  the  adsorption  complexes. 

Geometry  optimization  of  the  adsorption  complexes  involved  first  the  optimization  of  the 
surface  and  subsequently  optimization  of  the  adsorbed  sarin.  In  the  case  of  CaO,  inclusion  of  5 
atoms  on  the  surface  in  the  second  step  of  the  optimization  resulted  in  2.3  kcal/mol  increased 
binding  energy. 

Finally,  DFT/M062X  calculations  have  also  been  carried  out  for  sarin  on  MgO  (Mgi5Oi5) 
(cf.  (h)  in  Figure  3),  for  comparison  with  previous  work  on  this  system.6  For  these  calculations 
basis  sets  6-31G(d,p)  and  6-311G(d,p)  have  been  employed.  Addition  of  diffuse  functions  was 
not  possible  as  the  calculations  did  not  converge.  For  comparison,  the  previous  work  employed 
6-31G(d)  basis  set  with  DFT/B3LYP  and  MP2  calculations  on  different  models  of  sarin  on 

f: 

MgO.  All  the  DFT  and  TDDFT  calculations  of  the  present  work  have  been  carried  using 
Gaussian  09.28 


III.  Results  and  Discussions 
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3.1  Sarin  Geomerties 

Four  conformers  of  sarin  were  optimized  at  the  Cl  level.  The  minimum  energy  structure 
of  sarin  used  in  the  present  study  is  depicted  in  Figure  1.  This  is  the  geometry  labeled  sarin  II 
reported  in  Ref.  17.  The  energy  of  this  geometry  is  less  than  1  kcal/mol  above  the  minimum 
energy  geometry  reported  in  Ref.  17,  and  thus  for  the  purposes  of  this  research  is  considered 
unimportant.  The  remaining  conformers  differed  by  rotation  of  the  aliphatic  and  phosphate 
regions  and  changes  in  the  central  C-O-P  bond  angle. 

3.2  Binding  Geometries 

3.2.1  Overview 

In  Figure  3  the  structures  of  the  adsorption  complexes  obtained  by  DFT  and  Cl  are  given; 
in  Table  1  the  calculated  binding  energies  are  collected.  These  results  will  be  discussed  below 
starting  with  the  comparatively  simplest  case  of  sarin  on  CaO. 


(a)  Graphene,  DFT  (b)  Graphene,  DFT  (c)  Graphene,  Cl 
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(g)  CaO,  DFT  (h)  MgO,  DFT  (i)  CaO,  Cl 


Figure  3.  Optimum  geometries  of  the  adsorption  complexes  of  sarin  on  graphene,  (a)  and  (b) 
obtained  by  DFT  and  (c)  by  Cl,  on  graphane,  (d)  and  (e)  by  DFT  and  (f)  by  Cl,  on  CaO,  (g)  by 
DFT  and  (i)  by  Cl  and  on  MgO  (h)  by  DFT,  cf.  text. 
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3.2.2  Sarin  on  CaO 

Figures  3g-i  show  optimal  binding  geometries  of  sarin  to  oxide  surfaces.  The  binding 
appears  to  be  a  straight-forward  electrostatic  interaction.  The  P  atom  of  sarin  aligns  with  an  O 
atom  of  the  surface,  and  the  F  and  terminal  O  atoms  align  over  Ca  atoms,  creating  pairs  of 
opposite  charges.  The  center  of  mass  of  sarin  is  4.04  A  above  the  surface  in  the  Cl  calculations. 
In  all  calculations,  the  Ca-O/F  distance  is  on  the  order  of  2.4  A.  All  calculations  are  in  good 
agreement  on  the  binding  energy  (Table  1),  with  all  results  in  the  range  of  10-20  kcal/mol. 
Relaxation  of  the  CaO  lattice  during  binding  was  found  to  have  negligible  impact  on  the  binding 
energy  using  DFT.  For  comparison,  reported  values  of  binding  energy  of  sarin  on  MgO, 
calculated  by  B3LYP/6-31G(d)  and  a  large  model  (Mgi60i6)  is  2.9  kcal/mol  while  that  obtained 
by  MP2/6-31G(d)  and  a  small  cluster  (MgziOzO  is  50  kcal/mol.6  It  should  be  noted  that  in  the 
present  work  the  MgO  was  modelled  with  M062X  on  a  single  layer  of  Mgi50i5,  not  two  layers 
as  in  the  previous  work.  As  shown  in  Table  1,  significant  binding  of  20.0  kcal/mol  is  obtained 
for  sarin  on  MgO,  consistent  with  the  expectation  of  stronger  binding  of  sarin  to  MgO  than  CaO. 


Table  1.  Calculated  binding  energy  (BE)  of  sarin  at  different  surfaces  by  DFT  and  CL 


Surface  (structure  of  Figure  3)  /basis  set 

BE  (kcal/mol) 

BSSE-corrected 
BE  (kcal/mol) 

Graphene  (a)/  6-31G(d,p) 

6.7 

- 

Graphene  (b)/  6-31G(d,p) 

6.4 

2.6 

Graphene  (b)/  6-31 1+G(d,p) 

8.2 

4.8,  5.2* 

Graphane  (d)/  6-31G(d,p) 

9.0 

- 

Graphane  (e)/  6-31G(d,p) 

7.9 

4.9 

Graphane  (e)/  6-311+G(d,p) 

8.5 

6.4,  2.4* 

CaO  (g )/  6-31G(d,p) 

33.9 

16.1 

CaO  (g )/  6-311+G(d,p) 

22.3 

18.8,  13.2* 

MgO  (h)/6-31G(d,p) 

32.8 

19.3 

MgO  (h)/6-311G(d,p) 

31.8 

20.0 

Corresponding  Cl  value 


3.2.3  Sarin  on  graphene 

Unlike  the  results  for  CaO,  DFT  and  Cl  give  different  binding  geometries  for  sarin  on 
graphene  (Figure  3a-c).  Specifically,  DFT  models  predict  binding  through  the  sarin  oxygen  and 
fluorine  atoms,  while  Cl  predicts  a  flipped  orientation,  cf.  Figure  3b  vs  3c.  Consistency  between 
the  M062X/6-31G(d,p)  results  for  large  (Fig.  3a)  and  middle-sized  (Fig  3b)  graphene  surfaces 
suggests  that  the  M062X/6-31 1+G(d,p)  binding  energy  of  4.8  kcal/mol  should  be  taken  as 
accurate.  Despite  the  difference  in  orientation,  there  is  good  agreement  with  the  Cl  computed 
binding  energy  of  5.2  kcal/mol.  Cross-examinations  of  the  two  orientations  with  the  two  levels 
of  theory  indicate  that  the  energy  difference  between  the  two  orientations  is  small  (about  0. 1-0.4 
kcal/mol,  depending  on  the  model),  with  the  optimum  geometry  being  dependent  on  the 
particulars  of  the  method  used.  It  should  be  noted  that  in  both  orientations  the  dipole  moment  of 
sarin  is  directed  perpendicular  to  the  graphene  surface. 
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3.2.4  Sarin  on  graphane 

Of  the  three  surfaces  modelled,  graphane  exhibits  the  greatest  discrepancy  between  the 
DFT  and  Cl  results.  From  the  DFT  results,  sarin  binds  with  its  oxygen  atom  down  and  2.3  A 
from  a  nearest  surface  hydrogen  atom,  and  its  fluorine  atom  pointed  away  from  the  surface  (Fig. 
3d,e).  Predicted  binding  is  on  the  order  of  6  kcal/mol  -  slightly  greater  than  that  of  sarin  on 
graphene.  The  Cl  computations  suggest  a  binding  with  the  sarin  oxygen  and  fluorine  atoms  at 
the  surface  (Fig.  3f),  but  a  binding  of  only  2.4  kcal/mol.  Investigation  of  the  orientation  of  sarin 
on  graphane  obtained  by  Cl  employing  the  Cl  model  (C22H34),  as  well  as  the  larger  (C28H42) 
model,  using  the  DFT/M062X  method  show  that  for  the  smaller  Cl  model  only  the  Cl  structure 
(cf.  (f)  of  Figure  3)  is  obtained  whereas  with  the  larger  model  the  DFT  structure,  cf.  (e)  in  Figure 
3,  is  favored.  Therefore  there  is  a  significant  effect  of  the  size  of  the  model  employed  on  the 
calculated  orientation  of  the  adsorption  complex. 

3.3  Excited  States 

A  great  deal  of  experimental  and  theoretical  work  has  been  devoted  to  the  determination 
of  the  lowest  energy  structure  of  sarin  and  other  warfare  agents  in  the  ground  electronic  state,  but 
the  excited  states  have  not  been  studied  very  extensively.5  For  sarin,  a  low  intensity  absorption 
is  reported  with  the  threshold  electronic  excitation  occurring  in  the  vacuum-UV  (~7.1  eV  or 
higher)  with  varying  theoretical  estimates,  eg.  7.38  eV  obtained  by  TDDFT/B3LYP/6- 
311+G(2df,p)  and  9.9  eV  by  MP2/6-31+G(d)  calculations.9  Here  too  we  find  a  varying 
excitation  energy  for  the  lowest  singlet  excited  state  of  free  sarin,  depending  on  the  functional 
employed  in  TDDFT/6-311+G(d,p)  calculations  .  In  Figure  4  the  absorption  spectra  of  free  sarin 
obtained  by  TDDFT  calculations  and  ten  different  functionals,  offered  in  Gaussian  09,  are  shown 
The  lowest  onset  for  the  electronic  excitations  is  obtained  with  the  PBEPBE  functional  at  6.6  eV 
and  the  highest  with  M062X  at  8.1  eV.  Generally  low  oscillator  strengths  are  calculated  with  all 
functionals  employed.  Cl  calculations  on  the  excited  state  employing  the  basis  set  previously 
described  lead  to  an  excitation  energy  of  9.1  eV,  i.e.,  comparable  to  an  earlier  MP2/6-31+G(d) 
result.9  However,  inclusion  of  a  single  diffuse  s  function  (exponent  of  0.01)  in  the  phosphorus 
atom  basis  set  allows  an  accurate  Cl  estimate  of  the  lowest  Rydberg  excitations.  These  come  in 
at  8.0  eV  and  8.5  eV,  corresponding  to  excitations  from  the  oxygen/fluorine  lone-pair  manifold 
into  a  Rydberg  state.  This  is  in  good  agreement  with  the  present  CAM-B3LYP  and  M062X 
TDDFT  results  on  the  spectrum  of  free  sarin,  cf.  Figure  4. 

Next,  the  absorption  spectra  of  sarin  adsorbed  on  graphene  and  graphane  were  considered 
using  TDDFT/M062X  and  the  6-31 1+G(d,p)  basis  set.  For  the  graphene-sarin  system  the  lowest 
50  roots  were  calculated  in  the  region  1.98  eV  -  6.19  eV,  none  of  which  correspond  to  the  sarin- 
sarin  excitation.  The  lowest  transitions  correspond  to  graphene-graphene  excitations.  There  also 
exist  a  number  of  charge-transfer  type  excitations,  from  occupied  orbitals  of  graphene  to 
unoccupied  orbitals  of  sarin.  For  example,  states  calculated  at  3.37  and  4.2  eV  involve 
excitations  into  a  Rydberg  unoccupied  orbital  of  sarin.  Other  such  states  are  found  at  5.5,  6.08, 
6.18  and  6.19  eV.  The  calculations  on  graphane-sarin  involved  the  lowest  20  roots.  In  this  case 
the  lowest  unoccupied  orbital  or  LUMO  is  a  Rydberg  orbital  of  sarin  and  is  involved  in  several 
excited  states  of  the  complex.  In  particular,  the  lowest  excited  state  at  6.4  eV,  which  is 
characterized  by  the  HOMO— >LUMO  excitation,  corresponds  to  an  excitation  from  graphane  to 
sarin.  As  was  the  case  for  graphene-sarin,  the  calculations  did  not  reach  any  sarin-sarin  excited 
states.  It  should  be  noted  that  the  above  calculations  on  the  excited  states  of  the  adsorption 
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complexes,  just  as  in  the  case  of  free  sarin,  required  inclusion  of  diffuse  functions  in  the  basis  set. 
With  the  6-31G(d,p)  basis  set  the  first  25  roots  of  the  CaO-sarin  calculation  correspond  to 
excitations  involving  only  the  surface.  Calculations  with  the  6-311+G(d,p)  basis  on  the  CaO- 
sarin  complex  were  not  practical. 
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Figure  4.  Adsorption  spectra  of  free  sarin  calculated  by  TDDFT  employing  different 
functionals,  as  indicated. 


Cl  Computations  in  which  the  dominant  component  of  the  excitation  was  restricted  to  be 
into  the  Rydberg  orbital  of  sarin  were  used  to  determine  the  transition  energy  of  the  lowest  lying 
surface-to-sarin  electron  transfer  excitation.  For  graphane,  the  lowest  is  at  7.0  eV;  the  remainder 
are  above  the  free  sarin  excitation.  The  numerous  low-lying  virtual  orbitals  on  graphene  proved 
a  significant  challenge  for  isolating  the  charge-transfer  excitation,  and  the  best  estimate  places  it 
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near  5.3  eV.  Excitations  from  the  central  portion  of  the  CaO  model  surface  (excitations  from 
orbitals  localized  on  the  edges  of  the  model  were  disallowed)  provided  a  set  of  ten  electron- 
transfer  possibilities  ranging  in  transition  energy  from  6.2  eV  to  8.5  eV.  Therefore,  with  both 
methods,  the  existence  of  charge-transfer  excited  states  at  energies  lower  than  the  sarin-sarin 
excitation  is  indicated.  Similar  observations  have  been  reported  for  sarin  and  other 
organophosphorus  warfare  agents  adsorbed  on  a  Y-AI2O3  surface.5 

IV.  Summary 

Our  results  indicate  that  the  polarized  phosphate  region  of  the  sarin  molecule  plays  a 
critical  role  in  its  binding  to  aliphatic,  aromatic,  and  oxide  surfaces.  Specifically,  the  Lewis-base 
portion  acts  as  the  binding  agent.  For  the  surfaces  considered  in  the  present  study,  only  those 
that  are  ionic  (CaO)  or  have  large  polarizability  (graphene)  give  strong  binding  for  sarin. 

Binding  to  aliphatic  surfaces  is  weak.  Figure  5  depicts  Cl  potential  curves  for  interaction  of 
sarin  with  the  three  model  surfaces. 

E/kcalmol-1 


Figure  5.  Potential  curves  for  sarin  interaction  with  three  model  surfaces.  Infinite 
separation  limit  is  free  sarin  in  minumum  energy  geometry  and  bare  surface.  The  z-axis 
corresponds  to  the  height  of  the  sarin  center-of-ma 

ss  above  the  highest  atom  of  surface. 

It  is  concluded  that 

•  Sarin  binds  to  graphane,  graphene,  and  calcium  oxide  surfaces  by  2.4,  5.2,  and  13.2 
kcal/mol,  respectively  according  to  Cl,  with  corresponding  DFT  binding  energies  6.4,  4.8, 
and  18.8  kcal/mol,  respectively. 

•  Free  sarin  has  an  excitation  into  a  Rydberg  state  at  about  8  eV. 

•  Each  of  the  three  surfaces  studied  is  capable  of  transferring  an  electron  to  sarin  at  or 
below  the  sarin  Rydberg  excitation.  Both  DFT  and  many-electron  computational 
methods  gave  similar  results,  with  TDDFT  capable  of  determining  a  greater  total  number 
of  roots,  and  Cl  capable  to  giving  accurate  results  for  specific  states  of  interest. 
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